"bic.surv" <-
  
  function (x, ...) 
    
    UseMethod("bic.surv")

"bic.surv.data.frame" <-
  
  function (x, surv.t, cens, strict = FALSE, OR = 20, maxCol = 30, 
            
            prior.param = c(rep(0.5, ncol(x))), OR.fix = 2, nbest = 150, 
            
            factor.type = TRUE, factor.prior.adjust = FALSE, call = NULL, ...) 
    
  {
    
    leaps.bs <- function(info, coef, names.arg, nbest = nbest) {
      
      names.arg <- names.arg
      
      if (is.null(names.arg)) 
        
        names.arg <- c(as.character(1:9), LETTERS, letters)[1:ncol(info)]
      
      if (length(names.arg) < ncol(info)) 
        
        stop("Too few names")
      
      bIb <- coef %*% info %*% coef
      
      kx <- ncol(info)
      
      maxreg <- nbest * kx
      
      if (kx < 3) 
        
        stop("Too few independent variables")
      
      imeth <- 1
      
      df <- kx + 1
      
      Ib <- info %*% coef
      
      rr <- cbind(info, Ib)
      
      rr <- rbind(rr, c(Ib, bIb))
      
      it <- 0
      
      n.cols <- kx + 1
      
      nv <- kx + 1
      
      nf <- 0
      
      no <- 1e+05
      
      ib <- 1
      
      mb <- nbest
      
      nd <- n.cols
      
      nc <- 4 * n.cols
      
      rt <- matrix(rep(0, times = nd * nc), ncol = nc)
      
      rt[, 1:n.cols] <- rr
      
      iw <- c(1:(kx + 1), rep(0, times = 4 * nd))
      
      nw <- length(iw)
      
      rw <- rep(0, times = 2 * mb * kx + 7 * nd)
      
      nr <- length(rw)
      
      t1 <- 2
      
      s2 <- -1
      
      ne <- 0
      
      iv <- 0
      
      nret <- mb * kx
      
      Subss <- rep(0, times = nret)
      
      RSS <- Subss
      
      ans <- .Fortran("fwleaps", as.integer(nv), as.integer(it), 
                      
                      as.integer(kx), as.integer(nf), as.integer(no), as.integer(1), 
                      
                      as.double(2), as.integer(mb), as.double(rt), as.integer(nd), 
                      
                      as.integer(nc), as.integer(iw), as.integer(nw), as.double(rw), 
                      
                      as.integer(nr), as.double(t1), as.double(s2), as.integer(ne), 
                      
                      as.integer(iv), as.double(Subss), as.double(RSS), 
                      
                      as.integer(nret), PACKAGE = "BMA")
      
      regid <- ans[[21]]/2
      
      r2 <- ans[[20]]
      
      nreg <- sum(regid > 0)
      
      regid <- regid[1:nreg]
      
      r2 <- r2[1:nreg]
      
      which <- matrix(TRUE, nreg, kx)
      
      z <- regid
      
      which <- matrix(as.logical((rep.int(z, kx)%/%rep.int(2^((kx - 
                                                                 
                                                                 1):0), rep.int(length(z), kx)))%%2), byrow = FALSE, ncol = kx)
      
      size <- which %*% rep(1, kx)
      
      label <- character(nreg)
      
      sep <- if (all(nchar(names.arg) == 1)) 
        
        ""
      
      else ","
      
      for (i in 1:nreg) label[i] <- paste(names.arg[which[i, 
                                                          
                                                          ]], collapse = sep)
      
      ans <- list(r2 = r2, size = size, label = label, which = which)
      
      return(ans)
      
    }
    
    factor.names <- function(x) {
      
      out <- list()
      
      for (i in 1:ncol(x)) if (is.factor(x[, i])) 
        
        out[[i]] <- levels(x[, i])
      
      else out <- c(out, list(NULL))
      
      attributes(out)$names <- names(x)
      
      return(out)
      
    }
    
    dropcols <- function(x, surv.t, cens, maxCols = 30) {
      
      vnames <- attributes(x)$names
      
      nvar <- length(vnames)
      
      isfac <- rep(FALSE, times = nvar)
      
      for (i in 1:nvar) isfac[i] <- is.factor(x[, i])
      
      nlevels <- rep(NA, times = nvar)
      
      for (i in 1:nvar) if (isfac[i]) 
        
        nlevels[i] <- length(levels(x[, i]))
      
      any.dropped <- FALSE
      
      mm <- model.matrix(terms.formula(~., data = x), data = x)
      
      designx <- attributes(mm)$assign
      
      n.designx <- length(designx)
      
      designx.levels <- rep(1, times = n.designx)
      
      for (i in 2:n.designx) if (isfac[designx[i]]) 
        
        designx.levels[i] <- sum(designx[1:i] == designx[i]) + 
        
        1
      
      x.coxph <- data.frame(mm[, -1], surv.t = surv.t, cens = cens)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow", iter.max = 30)
      
      while (length(cox.out$coefficients) > maxCol) {
        
        any.dropped <- TRUE
        
        dropcox <- drop1(cox.out, test = "Chisq")
        
        dropped <- which.max(dropcox$"Pr(Chi)"[-1]) + 1
        
        x.coxph <- x.coxph[, -(dropped - 1)]
        
        designx.levels <- designx.levels[-dropped]
        
        designx <- designx[-dropped]
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         method = "breslow", iter.max = 30)
        
      }
      
      remaining.vars <- unique(designx[-1])
      
      new.nvar <- length(remaining.vars)
      
      dropped.vars <- vnames[-remaining.vars]
      
      dropped.levels <- NULL
      
      ncol.cox <- ncol(x.coxph) - 2
      
      x.coxph <- x.coxph[-((ncol.cox + 1):(ncol.cox + 2))]
      
      xx <- data.frame(matrix(rep(NA, times = new.nvar * nrow(x.coxph)), 
                              
                              ncol = new.nvar))
      
      new.names = rep(NA, times = new.nvar)
      
      for (i in 1:new.nvar) {
        
        cvar <- remaining.vars[i]
        
        lvls <- designx.levels[cvar == designx]
        
        if (isfac[cvar]) {
          
          if (length(lvls) != length(levels(x[, cvar]))) {
            
            newvar <- (as.matrix(x.coxph[, cvar == designx[-1]]) %*% 
                         
                         cbind(lvls - 1)) + 1
            
            xx[, i] <- factor(levels(x[, cvar])[newvar])
            
            new.names[i] <- vnames[cvar]
            
            removed.levels <- levels(x[, cvar])[-c(1, lvls)]
            
            dropped.levels <- c(dropped.levels, paste(vnames[cvar], 
                                                      
                                                      "_", removed.levels, sep = ""))
            
          }
          
          else {
            
            xx[, i] <- factor(x[, cvar])
            
            new.names[i] <- vnames[cvar]
            
          }
          
        }
        
        else {
          
          xx[, i] <- x[, cvar]
          
          new.names[i] <- vnames[cvar]
          
        }
        
      }
      
      dropped <- c(dropped.vars, dropped.levels)
      
      return(list(mm = xx, any.dropped = any.dropped, dropped = dropped, 
                  
                  var.names = new.names, remaining.vars = remaining.vars))
      
    }
    
    if (is.null(call)) 
      
      cl <- match.call()
    
    else cl <- call
    
    options(contrasts = c("contr.treatment", "contr.treatment"))
    
    x <- data.frame(x)
    
    prior.weight.denom <- 0.5^ncol(x)
    
    x.omit <- na.omit(x)
    
    used <- match(row.names(x), row.names(x.omit))
    
    omitted <- seq(nrow(x))[is.na(used)]
    
    if (length(omitted) > 0) {
      
      x <- x.omit
      
      surv.t <- surv.t[-omitted]
      
      cens <- cens[-omitted]
      
      warning(paste("There were ", length(omitted), "records deleted due to NA's"))
      
    }
    
    leaps.x <- x
    
    fac.levels <- rep(1, times = ncol(x))
    
    fn <- factor.names(x)
    
    factors <- !all(unlist(lapply(fn, is.null)))
    
    if (factors) {
      
      cdf <- cbind.data.frame(y = surv.t, x)
      
      mm <- model.matrix(formula(cdf), data = cdf)[, -1, drop = FALSE]
      
      mmm <- data.frame(matrix(mm, nrow = nrow(mm), byrow = FALSE))
      
      names(mmm) <- dimnames(mm)[[2]]
      
      output.names <- names(mmm)
      
      x.coxph <- data.frame(surv.t = surv.t, cens = cens, x)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow")
      
      cox.assign <- cox.out$assign
      
      for (tempi in length(cox.assign):1) {
        
        cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
        
        names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
        
      }
      
      cox.assign[[1]] <- numeric(0)
      
      names(cox.assign)[1] <- "(Intercept)"
      
      fac.levels <- unlist(lapply(cox.assign, length)[-1])
      
      if (factor.type) {
        
        for (i in 1:length(names(x))) {
          
          if (!is.null(fn[[i]])) {
            
            nx <- names(x)[i]
            
            coefs <- cox.out$coef[cox.out$assign[[i]]]
            
            old.vals <- x[, i]
            
            new.vals <- c(0, coefs)
            
            new.vec <- as.vector(new.vals[match(old.vals, 
                                                
                                                fn[[i]])])
            
            leaps.x[, nx] <- new.vec
            
          }
          
        }
        
      }
      
      else {
        
        fac.levels <- rep(1, times = ncol(mmm))
        
        new.prior <- NULL
        
        for (i in 1:length(names(x))) {
          
          addprior <- prior.param[i]
          
          if (!is.null(fn[[i]])) {
            
            k <- length(fn[[i]])
            
            if (factor.prior.adjust) 
              
              addprior <- rep(1 - (1 - prior.param[i])^(1/(k - 
                                                             
                                                             1)), k - 1)
            
            else addprior <- rep(prior.param[i], k - 1)
            
          }
          
          new.prior <- c(new.prior, addprior)
          
        }
        
        prior.param <- new.prior
        
        x <- leaps.x <- mmm
        
      }
      
    }
    
    
    
    
    
    xx <- data.frame()
    
    xx <- dropcols(leaps.x, surv.t, cens, maxCol)
    
    var.names <- xx$var.names
    
    remaining <- xx$remaining.vars
    
    leaps.x <- xx$mm
    
    reduced <- xx$any.dropped
    
    dropped <- NULL
    
    if (reduced) 
      
      dropped <- xx$dropped
    
    nvar <- length(x[1, ])
    
    x <- x[, remaining, drop = FALSE]
    
    x <- data.frame(x)
    
    fac.levels <- fac.levels[remaining]
    
    output.names <- list()
    
    for (i in 1:length(var.names)) {
      
      if (is.factor(x[, i])) 
        
        output.names[[i]] <- levels(x[, i])
      
      else output.names[[i]] <- NA
      
    }
    
    xnames <- names(x)
    
    names(leaps.x) <- var.names
    
    x.coxph <- data.frame(surv.t = surv.t, cens = cens, leaps.x)
    
    cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                     
                     method = "breslow")
    
    x.coxph.fac <- data.frame(surv.t = surv.t, cens = cens, x)
    
    cox.assign <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph.fac, 
                        
                        method = "breslow")$assign
    
    for (tempi in length(cox.assign):1) {
      
      cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
      
      names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
      
    }
    
    cox.assign[[1]] <- numeric(0)
    
    names(cox.assign)[1] <- "(Intercept)"
    
    n <- sum(cens)
    
    if (ncol(leaps.x) >= 3) {
      
      coef <- cox.out$coef
      
      info <- solve(cox.out$var)
      
      a <- leaps.bs(info, coef, names.arg = xnames, nbest = nbest)
      
      a$r2 <- c(0, a$r2)
      
      a$size <- c(0, a$size)
      
      a$label <- c("NULL", a$label)
      
      a$which <- rbind(rep(FALSE, ncol(x)), a$which)
      
      nmod <- length(a$size)
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(a$which * prior.mat + (!a$which) * (1 - 
                                                           
                                                           prior.mat), 1, prod)
      
      bIb <- as.numeric(coef %*% info %*% coef)
      
      lrt <- bIb - (a$r2 * bIb)
      
      bic <- lrt + (a$size) * log(n) - 2 * log(prior)
      
      occam <- bic - min(bic) < 2 * OR.fix * log(OR)
      
      size <- a$size[occam]
      
      label <- a$label[occam]
      
      which <- a$which[occam, , drop = FALSE]
      
      bic <- bic[occam]
      
      prior <- prior[occam]
      
    }
    
    else {
      
      nmod <- switch(ncol(x), 2, 4)
      
      bic <- label <- rep(0, nmod)
      
      which <- matrix(c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE), nmod, nmod/2)
      
      size <- c(0, 1, 1, 2)[1:nmod]
      
      sep <- ","
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(which * prior.mat + (!which) * (1 - prior.mat), 
                     
                     1, prod)
      
      for (k in 1:nmod) {
        
        if (k == 1) 
          
          label[k] <- "NULL"
        
        else label[k] <- paste(names(leaps.x)[which[k, ]], 
                               
                               collapse = sep)
        
      }
      
    }
    
    model.fits <- as.list(rep(0, length(label)))
    
    loglik <- rep(0, length(label))
    
    size <- rep(0, length(label))
    
    loglik.null <- cox.out$loglik[1]
    
    for (k in (1:length(label))) {
      
      if (sum(which[k, ]) != 0) {
        
        x.coxph <- data.frame(x[, which[k, ], drop = FALSE], 
                              
                              surv.t = surv.t, cens = cens)
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         iter.max = 30, method = "breslow")
        
        loglik[k] <- cox.out$loglik[2]
        
        size[k] <- length(cox.out$coef)
        
        model.fits[[k]] <- matrix(rep(0, 2 * length(cox.out$coef)), 
                                  
                                  ncol = 2)
        
        model.fits[[k]][, 1] <- cox.out$coef
        
        model.fits[[k]][, 2] <- sqrt(diag(cox.out$var))
        
      }
      
      else {
        
        loglik[k] <- loglik.null
        
      }
      
    }
    
    bic <- size * log(n) - 2 * (loglik - loglik.null) - 2 * log(prior)
    
    occam <- bic - min(bic) < 2 * log(OR)
    
    size <- size[occam]
    
    label <- label[occam]
    
    which <- which[occam, , drop = FALSE]
    
    bic <- bic[occam]
    
    prior <- prior[occam]
    
    model.fits <- model.fits[occam]
    
    postprob <- (exp(-0.5 * (bic - min(bic))))/sum(exp(-0.5 * 
                                                         
                                                         (bic - min(bic))))
    
    order.bic <- order(bic, size, label)
    
    size <- size[order.bic]
    
    label <- label[order.bic]
    
    which <- which[order.bic, , drop = FALSE]
    
    bic <- bic[order.bic]
    
    prior <- prior[order.bic]
    
    postprob <- postprob[order.bic]
    
    model.fits <- model.fits[order.bic]
    
    nmod <- length(size)
    
    if (strict & (nmod != 1)) {
      
      occam <- rep(TRUE, nmod)
      
      for (k in (2:nmod)) {
        
        for (j in (1:(k - 1))) {
          
          which.diff <- which[k, ] - which[j, ]
          
          if (all(which.diff >= 0)) 
            
            occam[k] <- FALSE
          
        }
        
      }
      
      size <- size[occam]
      
      label <- label[occam]
      
      which <- which[occam, , drop = FALSE]
      
      bic <- bic[occam]
      
      prior <- prior[occam]
      
      postprob <- postprob[occam]/(sum(postprob[occam]))
      
      model.fits <- model.fits[occam]
      
    }
    
    bic <- bic + 2 * log(prior)
    
    probne0 <- round(100 * t(which) %*% as.matrix(postprob), 
                     
                     1)
    
    nmod <- length(bic)
    
    nvar <- max(unlist(cox.assign))
    
    Ebi <- rep(0, nvar)
    
    SDbi <- rep(0, nvar)
    
    EbiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    sebiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    for (i in (1:ncol(x))) {
      
      whereisit <- cox.assign[[i + 1]]
      
      if (any(which[, i])) 
        
        for (k in (1:nmod)) if (which[k, i] == TRUE) {
          
          spot <- sum(which[k, (1:i)])
          
          posMk <- (c(0, cumsum(fac.levels[which[k, ]])) + 
                      
                      1)[spot]
          
          posMk <- posMk:(posMk + fac.levels[i] - 1)
          
          EbiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                 
                                                 1]
          
          sebiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                  
                                                  2]
          
        }
      
    }
    
    Ebi <- postprob %*% EbiMk
    
    Ebimat <- matrix(rep(Ebi, nmod), nrow = nmod, byrow = TRUE)
    
    SDbi <- sqrt(postprob %*% (sebiMk^2) + postprob %*% ((EbiMk - 
                                                            
                                                            Ebimat)^2))
    
    CSDbi <- rep(0, nvar)
    
    CEbi <- CSDbi
    
    for (i in (1:ncol(x))) {
      
      sel <- which[, i]
      
      if (sum(sel) > 0) {
        
        cpp <- rbind(postprob[sel]/sum(postprob[sel]))
        
        CEbi[cox.assign[[i + 1]]] <- as.numeric(cpp %*% EbiMk[sel, 
                                                              
                                                              cox.assign[[i + 1]]])
        
        CSDbi[cox.assign[[i + 1]]] <- sqrt(cpp %*% (sebiMk[sel, 
                                                           
                                                           cox.assign[[i + 1]]]^2) + cpp %*% ((EbiMk[sel, 
                                                                                                     
                                                                                                     cox.assign[[i + 1]]] - CEbi[cox.assign[[i + 1]]])^2))
        
      }
      
    }
    
    names(output.names) <- var.names
    
    result <- list(postprob = postprob, label = label, size = size, 
                   
                   bic = bic, prior.param = prior.param, prior.model.weights = prior/prior.weight.denom, 
                   
                   which = which, probne0 = c(probne0), postmean = as.vector(Ebi), 
                   
                   postsd = as.vector(SDbi), condpostmean = CEbi, condpostsd = CSDbi, 
                   
                   mle = EbiMk, se = sebiMk, namesx = var.names, reduced = reduced, 
                   
                   dropped = dropped, call = cl, n.models = length(postprob), 
                   
                   n.vars = length(probne0), nests = length(Ebi), output.names = output.names, 
                   
                   assign = cox.assign, factor.type = factor.type)
    
    class(result) <- "bic.surv"
    
    result
    
  }

"bic.surv.formula" <-
  
  function (f, data, strict = FALSE, OR = 20, maxCol = 30, prior.param = c(rep(0.5, 
                                                                               
                                                                               ncol(x))), OR.fix = 2, nbest = 150, factor.type = TRUE, factor.prior.adjust = FALSE, 
            
            call = NULL, ...) 
    
  {
    
    cl <- match.call()
    
    tms <- terms(f, data = data)
    
    fmatrix <- attr(tms, "factors")
    
    tms.order <- attr(tms, "order")
    
    tms.labels <- attr(tms, "term.labels")
    
    mm <- model.matrix(tms, data = data)
    
    assn <- attr(mm, "assign")
    
    nterms <- max(assn)
    
    datalist <- eval(attr(tms, "variables"), envir = data)
    
    nvar <- nrow(fmatrix) - 1
    
    isvarfac <- rep(NA, times = nvar)
    
    for (i in 1:nvar) isvarfac[i] <- is.factor(datalist[[i + 
                                                           
                                                           1]])
    
    istermfac <- rep(NA, times = nterms)
    
    for (i in 1:nterms) {
      
      cterms <- fmatrix[-1, i] == 1
      
      istermfac[i] <- sum(isvarfac[cterms] == FALSE) == 0
      
    }
    
    surv.t.name <- all.vars(f)[1]
    
    cens.name <- all.vars(f)[2]
    
    moddata <- data.frame(rep(NA, times = dim(mm)[1]))
    
    cnames <- NULL
    
    for (i in 1:nterms) {
      
      if (istermfac[i]) {
        
        if (tms.order[i] == 1) {
          
          moddata <- cbind(moddata, datalist[[i + 1]])
          
          cnames <- c(cnames, tms.labels[i])
          
        }
        
        else {
          
          sel <- assn == i
          
          nlev <- sum(sel)
          
          newfac.index <- (mm[, sel] %*% cbind(1:nlev)) + 
            
            1
          
          facnames <- c("ref", colnames(mm)[sel])
          
          newfac <- facnames[newfac.index]
          
          newfac <- factor(newfac)
          
          moddata <- cbind(moddata, newfac)
          
          cnames <- c(cnames, paste(tms.labels[i], "..", 
                                    
                                    sep = ""))
          
        }
        
      }
      
      else {
        
        sel <- assn == i
        
        moddata <- cbind(moddata, mm[, sel])
        
        cnames <- c(cnames, colnames(mm)[sel])
        
      }
      
    }
    
    moddata <- moddata[, -1]
    
    cnames <- gsub(":", ".", cnames)
    
    moddata <- cbind(moddata, datalist[[1]][, 1], datalist[[1]][, 
                                                                
                                                                2])
    
    colnames(moddata) <- c(cnames, surv.t.name, cens.name)
    
    nv <- ncol(moddata) - 2
    
    surv.t <- moddata[, nv + 1]
    
    cens <- moddata[, nv + 2]
    
    x <- moddata[, 1:nv]
    
    bic.surv(x, surv.t, cens, strict = FALSE, OR = OR, maxCol = maxCol, 
             
             prior.param = prior.param, OR.fix = OR.fix, nbest = nbest, 
             
             factor.type = factor.type, factor.prior.adjust = factor.prior.adjust, 
             
             call = cl)
    
  }

"bic.surv.matrix" <-
  
  function (x, surv.t, cens, strict = FALSE, OR = 20, maxCol = 30, 
            
            prior.param = c(rep(0.5, ncol(x))), OR.fix = 2, nbest = 150, 
            
            factor.type = TRUE, factor.prior.adjust = FALSE, call = NULL, ...) 
    
  {
    
    leaps.bs <- function(info, coef, names.arg, nbest = nbest) {
      
      names.arg <- names.arg
      
      if (is.null(names.arg)) 
        
        names.arg <- c(as.character(1:9), LETTERS, letters)[1:ncol(info)]
      
      if (length(names.arg) < ncol(info)) 
        
        stop("Too few names")
      
      bIb <- coef %*% info %*% coef
      
      kx <- ncol(info)
      
      maxreg <- nbest * kx
      
      if (kx < 3) 
        
        stop("Too few independent variables")
      
      imeth <- 1
      
      df <- kx + 1
      
      Ib <- info %*% coef
      
      rr <- cbind(info, Ib)
      
      rr <- rbind(rr, c(Ib, bIb))
      
      it <- 0
      
      n.cols <- kx + 1
      
      nv <- kx + 1
      
      nf <- 0
      
      no <- 1e+05
      
      ib <- 1
      
      mb <- nbest
      
      nd <- n.cols
      
      nc <- 4 * n.cols
      
      rt <- matrix(rep(0, times = nd * nc), ncol = nc)
      
      rt[, 1:n.cols] <- rr
      
      iw <- c(1:(kx + 1), rep(0, times = 4 * nd))
      
      nw <- length(iw)
      
      rw <- rep(0, times = 2 * mb * kx + 7 * nd)
      
      nr <- length(rw)
      
      t1 <- 2
      
      s2 <- -1
      
      ne <- 0
      
      iv <- 0
      
      nret <- mb * kx
      
      Subss <- rep(0, times = nret)
      
      RSS <- Subss
      
      ans <- .Fortran("fwleaps", as.integer(nv), as.integer(it), 
                      
                      as.integer(kx), as.integer(nf), as.integer(no), as.integer(1), 
                      
                      as.double(2), as.integer(mb), as.double(rt), as.integer(nd), 
                      
                      as.integer(nc), as.integer(iw), as.integer(nw), as.double(rw), 
                      
                      as.integer(nr), as.double(t1), as.double(s2), as.integer(ne), 
                      
                      as.integer(iv), as.double(Subss), as.double(RSS), 
                      
                      as.integer(nret), PACKAGE = "BMA")
      
      regid <- ans[[21]]/2
      
      r2 <- ans[[20]]
      
      nreg <- sum(regid > 0)
      
      regid <- regid[1:nreg]
      
      r2 <- r2[1:nreg]
      
      which <- matrix(TRUE, nreg, kx)
      
      z <- regid
      
      which <- matrix(as.logical((rep.int(z, kx)%/%rep.int(2^((kx - 
                                                                 
                                                                 1):0), rep.int(length(z), kx)))%%2), byrow = FALSE, ncol = kx)
      
      size <- which %*% rep(1, kx)
      
      label <- character(nreg)
      
      sep <- if (all(nchar(names.arg) == 1)) 
        
        ""
      
      else ","
      
      for (i in 1:nreg) label[i] <- paste(names.arg[which[i, 
                                                          
                                                          ]], collapse = sep)
      
      ans <- list(r2 = r2, size = size, label = label, which = which)
      
      return(ans)
      
    }
    
    factor.names <- function(x) {
      
      out <- list()
      
      for (i in 1:ncol(x)) if (is.factor(x[, i])) 
        
        out[[i]] <- levels(x[, i])
      
      else out <- c(out, list(NULL))
      
      attributes(out)$names <- names(x)
      
      return(out)
      
    }
    
    dropcols <- function(x, surv.t, cens, maxCols = 30) {
      
      vnames <- attributes(x)$names
      
      nvar <- length(vnames)
      
      isfac <- rep(FALSE, times = nvar)
      
      for (i in 1:nvar) isfac[i] <- is.factor(x[, i])
      
      nlevels <- rep(NA, times = nvar)
      
      for (i in 1:nvar) if (isfac[i]) 
        
        nlevels[i] <- length(levels(x[, i]))
      
      any.dropped <- FALSE
      
      mm <- model.matrix(terms.formula(~., data = x), data = x)
      
      designx <- attributes(mm)$assign
      
      n.designx <- length(designx)
      
      designx.levels <- rep(1, times = n.designx)
      
      for (i in 2:n.designx) if (isfac[designx[i]]) 
        
        designx.levels[i] <- sum(designx[1:i] == designx[i]) + 
        
        1
      
      x.coxph <- data.frame(mm[, -1], surv.t = surv.t, cens = cens)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow", iter.max = 30)
      
      while (length(cox.out$coefficients) > maxCol) {
        
        any.dropped <- TRUE
        
        dropcox <- drop1(cox.out, test = "Chisq")
        
        dropped <- which.max(dropcox$"Pr(Chi)"[-1]) + 1
        
        x.coxph <- x.coxph[, -(dropped - 1)]
        
        designx.levels <- designx.levels[-dropped]
        
        designx <- designx[-dropped]
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         method = "breslow", iter.max = 30)
        
      }
      
      remaining.vars <- unique(designx[-1])
      
      new.nvar <- length(remaining.vars)
      
      dropped.vars <- vnames[-remaining.vars]
      
      dropped.levels <- NULL
      
      ncol.cox <- ncol(x.coxph) - 2
      
      x.coxph <- x.coxph[-((ncol.cox + 1):(ncol.cox + 2))]
      
      xx <- data.frame(matrix(rep(NA, times = new.nvar * nrow(x.coxph)), 
                              
                              ncol = new.nvar))
      
      new.names = rep(NA, times = new.nvar)
      
      for (i in 1:new.nvar) {
        
        cvar <- remaining.vars[i]
        
        lvls <- designx.levels[cvar == designx]
        
        if (isfac[cvar]) {
          
          if (length(lvls) != length(levels(x[, cvar]))) {
            
            newvar <- (as.matrix(x.coxph[, cvar == designx[-1]]) %*% 
                         
                         cbind(lvls - 1)) + 1
            
            xx[, i] <- factor(levels(x[, cvar])[newvar])
            
            new.names[i] <- vnames[cvar]
            
            removed.levels <- levels(x[, cvar])[-c(1, lvls)]
            
            dropped.levels <- c(dropped.levels, paste(vnames[cvar], 
                                                      
                                                      "_", removed.levels, sep = ""))
            
          }
          
          else {
            
            xx[, i] <- factor(x[, cvar])
            
            new.names[i] <- vnames[cvar]
            
          }
          
        }
        
        else {
          
          xx[, i] <- x[, cvar]
          
          new.names[i] <- vnames[cvar]
          
        }
        
      }
      
      dropped <- c(dropped.vars, dropped.levels)
      
      return(list(mm = xx, any.dropped = any.dropped, dropped = dropped, 
                  
                  var.names = new.names, remaining.vars = remaining.vars))
      
    }
    
    if (is.null(call)) 
      
      cl <- match.call()
    
    else cl <- call
    
    options(contrasts = c("contr.treatment", "contr.treatment"))
    
    x <- data.frame(x)
    
    prior.weight.denom <- 0.5^ncol(x)
    
    x.omit <- na.omit(x)
    
    used <- match(row.names(x), row.names(x.omit))
    
    omitted <- seq(nrow(x))[is.na(used)]
    
    if (length(omitted) > 0) {
      
      x <- x.omit
      
      surv.t <- surv.t[-omitted]
      
      cens <- cens[-omitted]
      
      warning(paste("There were ", length(omitted), "records deleted due to NA's"))
      
    }
    
    leaps.x <- x
    
    fac.levels <- rep(1, times = ncol(x))
    
    fn <- factor.names(x)
    
    factors <- !all(unlist(lapply(fn, is.null)))
    
    if (factors) {
      
      cdf <- cbind.data.frame(y = surv.t, x)
      
      mm <- model.matrix(formula(cdf), data = cdf)[, -1, drop = FALSE]
      
      mmm <- data.frame(matrix(mm, nrow = nrow(mm), byrow = FALSE))
      
      names(mmm) <- dimnames(mm)[[2]]
      
      output.names <- names(mmm)
      
      x.coxph <- data.frame(surv.t = surv.t, cens = cens, x)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow")
      
      cox.assign <- cox.out$assign
      
      for (tempi in length(cox.assign):1) {
        
        cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
        
        names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
        
      }
      
      cox.assign[[1]] <- numeric(0)
      
      names(cox.assign)[1] <- "(Intercept)"
      
      fac.levels <- unlist(lapply(cox.assign, length)[-1])
      
      if (factor.type) {
        
        for (i in 1:length(names(x))) {
          
          if (!is.null(fn[[i]])) {
            
            nx <- names(x)[i]
            
            coefs <- cox.out$coef[cox.out$assign[[i]]]
            
            old.vals <- x[, i]
            
            new.vals <- c(0, coefs)
            
            new.vec <- as.vector(new.vals[match(old.vals, 
                                                
                                                fn[[i]])])
            
            leaps.x[, nx] <- new.vec
            
          }
          
        }
        
      }
      
      else {
        
        fac.levels <- rep(1, times = ncol(mmm))
        
        new.prior <- NULL
        
        for (i in 1:length(names(x))) {
          
          addprior <- prior.param[i]
          
          if (!is.null(fn[[i]])) {
            
            k <- length(fn[[i]])
            
            if (factor.prior.adjust) 
              
              addprior <- rep(1 - (1 - prior.param[i])^(1/(k - 
                                                             
                                                             1)), k - 1)
            
            else addprior <- rep(prior.param[i], k - 1)
            
          }
          
          new.prior <- c(new.prior, addprior)
          
        }
        
        prior.param <- new.prior
        
        x <- leaps.x <- mmm
        
      }
      
    }
    
    xx <- data.frame()
    
    xx <- dropcols(leaps.x, surv.t, cens, maxCol)
    
    var.names <- xx$var.names
    
    remaining <- xx$remaining.vars
    
    leaps.x <- xx$mm
    
    reduced <- xx$any.dropped
    
    dropped <- NULL
    
    if (reduced) 
      
      dropped <- xx$dropped
    
    nvar <- length(x[1, ])
    
    x <- x[, remaining, drop = FALSE]
    
    x <- data.frame(x)
    
    fac.levels <- fac.levels[remaining]
    
    output.names <- list()
    
    for (i in 1:length(var.names)) {
      
      if (is.factor(x[, i])) 
        
        output.names[[i]] <- levels(x[, i])
      
      else output.names[[i]] <- NA
      
    }
    
    xnames <- names(x)
    
    names(leaps.x) <- var.names
    
    x.coxph <- data.frame(surv.t = surv.t, cens = cens, leaps.x)
    
    cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                     
                     method = "breslow")
    
    x.coxph.fac <- data.frame(surv.t = surv.t, cens = cens, x)
    
    cox.assign <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph.fac, 
                        
                        method = "breslow")$assign
    
    for (tempi in length(cox.assign):1) {
      
      cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
      
      names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
      
    }
    
    cox.assign[[1]] <- numeric(0)
    
    names(cox.assign)[1] <- "(Intercept)"
    
    n <- sum(cens)
    
    if (ncol(leaps.x) >= 3) {
      
      coef <- cox.out$coef
      
      info <- solve(cox.out$var)
      
      a <- leaps.bs(info, coef, names.arg = xnames, nbest = nbest)
      
      a$r2 <- c(0, a$r2)
      
      a$size <- c(0, a$size)
      
      a$label <- c("NULL", a$label)
      
      a$which <- rbind(rep(FALSE, ncol(x)), a$which)
      
      nmod <- length(a$size)
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(a$which * prior.mat + (!a$which) * (1 - 
                                                           
                                                           prior.mat), 1, prod)
      
      bIb <- as.numeric(coef %*% info %*% coef)
      
      lrt <- bIb - (a$r2 * bIb)
      
      bic <- lrt + (a$size) * log(n) - 2 * log(prior)
      
      occam <- bic - min(bic) < 2 * OR.fix * log(OR)
      
      size <- a$size[occam]
      
      label <- a$label[occam]
      
      which <- a$which[occam, , drop = FALSE]
      
      bic <- bic[occam]
      
      prior <- prior[occam]
      
    }
    
    else {
      
      nmod <- switch(ncol(x), 2, 4)
      
      bic <- label <- rep(0, nmod)
      
      which <- matrix(c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE), nmod, nmod/2)
      
      size <- c(0, 1, 1, 2)[1:nmod]
      
      sep <- ","
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(which * prior.mat + (!which) * (1 - prior.mat), 
                     
                     1, prod)
      
      for (k in 1:nmod) {
        
        if (k == 1) 
          
          label[k] <- "NULL"
        
        else label[k] <- paste(names(leaps.x)[which[k, ]], 
                               
                               collapse = sep)
        
      }
      
    }
    
    model.fits <- as.list(rep(0, length(label)))
    
    loglik <- rep(0, length(label))
    
    size <- rep(0, length(label))
    
    loglik.null <- cox.out$loglik[1]
    
    for (k in (1:length(label))) {
      
      if (sum(which[k, ]) != 0) {
        
        x.coxph <- data.frame(x[, which[k, ], drop = FALSE], 
                              
                              surv.t = surv.t, cens = cens)
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         iter.max = 30, method = "breslow")
        
        loglik[k] <- cox.out$loglik[2]
        
        size[k] <- length(cox.out$coef)
        
        model.fits[[k]] <- matrix(rep(0, 2 * length(cox.out$coef)), 
                                  
                                  ncol = 2)
        
        model.fits[[k]][, 1] <- cox.out$coef
        
        model.fits[[k]][, 2] <- sqrt(diag(cox.out$var))
        
      }
      
      else {
        
        loglik[k] <- loglik.null
        
      }
      
    }
    
    bic <- size * log(n) - 2 * (loglik - loglik.null) - 2 * log(prior)
    
    occam <- bic - min(bic) < 2 * log(OR)
    
    size <- size[occam]
    
    label <- label[occam]
    
    which <- which[occam, , drop = FALSE]
    
    bic <- bic[occam]
    
    prior <- prior[occam]
    
    model.fits <- model.fits[occam]
    
    postprob <- (exp(-0.5 * (bic - min(bic))))/sum(exp(-0.5 * 
                                                         
                                                         (bic - min(bic))))
    
    order.bic <- order(bic, size, label)
    
    size <- size[order.bic]
    
    label <- label[order.bic]
    
    which <- which[order.bic, , drop = FALSE]
    
    bic <- bic[order.bic]
    
    prior <- prior[order.bic]
    
    postprob <- postprob[order.bic]
    
    model.fits <- model.fits[order.bic]
    
    nmod <- length(size)
    
    if (strict & (nmod != 1)) {
      
      occam <- rep(TRUE, nmod)
      
      for (k in (2:nmod)) {
        
        for (j in (1:(k - 1))) {
          
          which.diff <- which[k, ] - which[j, ]
          
          if (all(which.diff >= 0)) 
            
            occam[k] <- FALSE
          
        }
        
      }
      
      size <- size[occam]
      
      label <- label[occam]
      
      which <- which[occam, , drop = FALSE]
      
      bic <- bic[occam]
      
      prior <- prior[occam]
      
      postprob <- postprob[occam]/(sum(postprob[occam]))
      
      model.fits <- model.fits[occam]
      
    }
    
    bic <- bic + 2 * log(prior)
    
    probne0 <- round(100 * t(which) %*% as.matrix(postprob), 
                     
                     1)
    
    nmod <- length(bic)
    
    nvar <- max(unlist(cox.assign))
    
    Ebi <- rep(0, nvar)
    
    SDbi <- rep(0, nvar)
    
    EbiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    sebiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    for (i in (1:ncol(x))) {
      
      whereisit <- cox.assign[[i + 1]]
      
      if (any(which[, i])) 
        
        for (k in (1:nmod)) if (which[k, i] == TRUE) {
          
          spot <- sum(which[k, (1:i)])
          
          posMk <- (c(0, cumsum(fac.levels[which[k, ]])) + 
                      
                      1)[spot]
          
          posMk <- posMk:(posMk + fac.levels[i] - 1)
          
          EbiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                 
                                                 1]
          
          sebiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                  
                                                  2]
          
        }
      
    }
    
    Ebi <- postprob %*% EbiMk
    
    Ebimat <- matrix(rep(Ebi, nmod), nrow = nmod, byrow = TRUE)
    
    SDbi <- sqrt(postprob %*% (sebiMk^2) + postprob %*% ((EbiMk - 
                                                            
                                                            Ebimat)^2))
    
    CSDbi <- rep(0, nvar)
    
    CEbi <- CSDbi
    
    for (i in (1:ncol(x))) {
      
      sel <- which[, i]
      
      if (sum(sel) > 0) {
        
        cpp <- rbind(postprob[sel]/sum(postprob[sel]))
        
        CEbi[cox.assign[[i + 1]]] <- as.numeric(cpp %*% EbiMk[sel, 
                                                              
                                                              cox.assign[[i + 1]]])
        
        CSDbi[cox.assign[[i + 1]]] <- sqrt(cpp %*% (sebiMk[sel, 
                                                           
                                                           cox.assign[[i + 1]]]^2) + cpp %*% ((EbiMk[sel, 
                                                                                                     
                                                                                                     cox.assign[[i + 1]]] - CEbi[cox.assign[[i + 1]]])^2))
        
      }
      
    }
    
    names(output.names) <- var.names
    
    result <- list(postprob = postprob, label = label, size = size, 
                   
                   bic = bic, prior.param = prior.param, prior.model.weights = prior/prior.weight.denom, 
                   
                   which = which, probne0 = c(probne0), postmean = as.vector(Ebi), 
                   
                   postsd = as.vector(SDbi), condpostmean = CEbi, condpostsd = CSDbi, 
                   
                   mle = EbiMk, se = sebiMk, namesx = var.names, reduced = reduced, 
                   
                   dropped = dropped, call = cl, n.models = length(postprob), 
                   
                   n.vars = length(probne0), nests = length(Ebi), output.names = output.names, 
                   
                   assign = cox.assign, factor.type = factor.type)
    
    class(result) <- "bic.surv"
    
    result
    
  }


"aic.surv" <-
  function (x, ...) 
    UseMethod("aic.surv")

"aic.surv.data.frame" <-
  
  function (x, surv.t, cens, strict = FALSE, OR = 20, maxCol = 30, 
            
            prior.param = c(rep(0.5, ncol(x))), OR.fix = 2, nbest = 150, 
            
            factor.type = TRUE, factor.prior.adjust = FALSE, call = NULL, ...) 
    
  {
    
    leaps.bs <- function(info, coef, names.arg, nbest = nbest) {
      
      names.arg <- names.arg
      
      if (is.null(names.arg)) 
        
        names.arg <- c(as.character(1:9), LETTERS, letters)[1:ncol(info)]
      
      if (length(names.arg) < ncol(info)) 
        
        stop("Too few names")
      
      bIb <- coef %*% info %*% coef
      
      kx <- ncol(info)
      
      maxreg <- nbest * kx
      
      if (kx < 3) 
        
        stop("Too few independent variables")
      
      imeth <- 1
      
      df <- kx + 1
      
      Ib <- info %*% coef
      
      rr <- cbind(info, Ib)
      
      rr <- rbind(rr, c(Ib, bIb))
      
      it <- 0
      
      n.cols <- kx + 1
      
      nv <- kx + 1
      
      nf <- 0
      
      no <- 1e+05
      
      ib <- 1
      
      mb <- nbest
      
      nd <- n.cols
      
      nc <- 4 * n.cols
      
      rt <- matrix(rep(0, times = nd * nc), ncol = nc)
      
      rt[, 1:n.cols] <- rr
      
      iw <- c(1:(kx + 1), rep(0, times = 4 * nd))
      
      nw <- length(iw)
      
      rw <- rep(0, times = 2 * mb * kx + 7 * nd)
      
      nr <- length(rw)
      
      t1 <- 2
      
      s2 <- -1
      
      ne <- 0
      
      iv <- 0
      
      nret <- mb * kx
      
      Subss <- rep(0, times = nret)
      
      RSS <- Subss
      
      ans <- .Fortran("fwleaps", as.integer(nv), as.integer(it), 
                      
                      as.integer(kx), as.integer(nf), as.integer(no), as.integer(1), 
                      
                      as.double(2), as.integer(mb), as.double(rt), as.integer(nd), 
                      
                      as.integer(nc), as.integer(iw), as.integer(nw), as.double(rw), 
                      
                      as.integer(nr), as.double(t1), as.double(s2), as.integer(ne), 
                      
                      as.integer(iv), as.double(Subss), as.double(RSS), 
                      
                      as.integer(nret), PACKAGE = "BMA")
      
      regid <- ans[[21]]/2
      
      r2 <- ans[[20]]
      
      nreg <- sum(regid > 0)
      
      regid <- regid[1:nreg]
      
      r2 <- r2[1:nreg]
      
      which <- matrix(TRUE, nreg, kx)
      
      z <- regid
      
      which <- matrix(as.logical((rep.int(z, kx)%/%rep.int(2^((kx - 
                                                                 
                                                                 1):0), rep.int(length(z), kx)))%%2), byrow = FALSE, ncol = kx)
      
      size <- which %*% rep(1, kx)
      
      label <- character(nreg)
      
      sep <- if (all(nchar(names.arg) == 1)) 
        
        ""
      
      else ","
      
      for (i in 1:nreg) label[i] <- paste(names.arg[which[i, 
                                                          
                                                          ]], collapse = sep)
      
      ans <- list(r2 = r2, size = size, label = label, which = which)
      
      return(ans)
      
    }
    
    factor.names <- function(x) {
      
      out <- list()
      
      for (i in 1:ncol(x)) if (is.factor(x[, i])) 
        
        out[[i]] <- levels(x[, i])
      
      else out <- c(out, list(NULL))
      
      attributes(out)$names <- names(x)
      
      return(out)
      
    }
    
    dropcols <- function(x, surv.t, cens, maxCols = 30) {
      
      vnames <- attributes(x)$names
      
      nvar <- length(vnames)
      
      isfac <- rep(FALSE, times = nvar)
      
      for (i in 1:nvar) isfac[i] <- is.factor(x[, i])
      
      nlevels <- rep(NA, times = nvar)
      
      for (i in 1:nvar) if (isfac[i]) 
        
        nlevels[i] <- length(levels(x[, i]))
      
      any.dropped <- FALSE
      
      mm <- model.matrix(terms.formula(~., data = x), data = x)
      
      designx <- attributes(mm)$assign
      
      n.designx <- length(designx)
      
      designx.levels <- rep(1, times = n.designx)
      
      for (i in 2:n.designx) if (isfac[designx[i]]) 
        
        designx.levels[i] <- sum(designx[1:i] == designx[i]) + 
        
        1
      
      x.coxph <- data.frame(mm[, -1], surv.t = surv.t, cens = cens)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow", iter.max = 30)
      
      while (length(cox.out$coefficients) > maxCol) {
        
        any.dropped <- TRUE
        
        dropcox <- drop1(cox.out, test = "Chisq")
        
        dropped <- which.max(dropcox$"Pr(Chi)"[-1]) + 1
        
        x.coxph <- x.coxph[, -(dropped - 1)]
        
        designx.levels <- designx.levels[-dropped]
        
        designx <- designx[-dropped]
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         method = "breslow", iter.max = 30)
        
      }
      
      remaining.vars <- unique(designx[-1])
      
      new.nvar <- length(remaining.vars)
      
      dropped.vars <- vnames[-remaining.vars]
      
      dropped.levels <- NULL
      
      ncol.cox <- ncol(x.coxph) - 2
      
      x.coxph <- x.coxph[-((ncol.cox + 1):(ncol.cox + 2))]
      
      xx <- data.frame(matrix(rep(NA, times = new.nvar * nrow(x.coxph)), 
                              
                              ncol = new.nvar))
      
      new.names = rep(NA, times = new.nvar)
      
      for (i in 1:new.nvar) {
        
        cvar <- remaining.vars[i]
        
        lvls <- designx.levels[cvar == designx]
        
        if (isfac[cvar]) {
          
          if (length(lvls) != length(levels(x[, cvar]))) {
            
            newvar <- (as.matrix(x.coxph[, cvar == designx[-1]]) %*% 
                         
                         cbind(lvls - 1)) + 1
            
            xx[, i] <- factor(levels(x[, cvar])[newvar])
            
            new.names[i] <- vnames[cvar]
            
            removed.levels <- levels(x[, cvar])[-c(1, lvls)]
            
            dropped.levels <- c(dropped.levels, paste(vnames[cvar], 
                                                      
                                                      "_", removed.levels, sep = ""))
            
          }
          
          else {
            
            xx[, i] <- factor(x[, cvar])
            
            new.names[i] <- vnames[cvar]
            
          }
          
        }
        
        else {
          
          xx[, i] <- x[, cvar]
          
          new.names[i] <- vnames[cvar]
          
        }
        
      }
      
      dropped <- c(dropped.vars, dropped.levels)
      
      return(list(mm = xx, any.dropped = any.dropped, dropped = dropped, 
                  
                  var.names = new.names, remaining.vars = remaining.vars))
      
    }
    
    if (is.null(call)) 
      
      cl <- match.call()
    
    else cl <- call
    
    options(contrasts = c("contr.treatment", "contr.treatment"))
    
    x <- data.frame(x)
    
    prior.weight.denom <- 0.5^ncol(x)
    
    x.omit <- na.omit(x)
    
    used <- match(row.names(x), row.names(x.omit))
    
    omitted <- seq(nrow(x))[is.na(used)]
    
    if (length(omitted) > 0) {
      
      x <- x.omit
      
      surv.t <- surv.t[-omitted]
      
      cens <- cens[-omitted]
      
      warning(paste("There were ", length(omitted), "records deleted due to NA's"))
      
    }
    
    leaps.x <- x
    
    fac.levels <- rep(1, times = ncol(x))
    
    fn <- factor.names(x)
    
    factors <- !all(unlist(lapply(fn, is.null)))
    
    if (factors) {
      
      cdf <- cbind.data.frame(y = surv.t, x)
      
      mm <- model.matrix(formula(cdf), data = cdf)[, -1, drop = FALSE]
      
      mmm <- data.frame(matrix(mm, nrow = nrow(mm), byrow = FALSE))
      
      names(mmm) <- dimnames(mm)[[2]]
      
      output.names <- names(mmm)
      
      x.coxph <- data.frame(surv.t = surv.t, cens = cens, x)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow")
      
      cox.assign <- cox.out$assign
      
      for (tempi in length(cox.assign):1) {
        
        cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
        
        names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
        
      }
      
      cox.assign[[1]] <- numeric(0)
      
      names(cox.assign)[1] <- "(Intercept)"
      
      fac.levels <- unlist(lapply(cox.assign, length)[-1])
      
      if (factor.type) {
        
        for (i in 1:length(names(x))) {
          
          if (!is.null(fn[[i]])) {
            
            nx <- names(x)[i]
            
            coefs <- cox.out$coef[cox.out$assign[[i]]]
            
            old.vals <- x[, i]
            
            new.vals <- c(0, coefs)
            
            new.vec <- as.vector(new.vals[match(old.vals, 
                                                
                                                fn[[i]])])
            
            leaps.x[, nx] <- new.vec
            
          }
          
        }
        
      }
      
      else {
        
        fac.levels <- rep(1, times = ncol(mmm))
        
        new.prior <- NULL
        
        for (i in 1:length(names(x))) {
          
          addprior <- prior.param[i]
          
          if (!is.null(fn[[i]])) {
            
            k <- length(fn[[i]])
            
            if (factor.prior.adjust) 
              
              addprior <- rep(1 - (1 - prior.param[i])^(1/(k - 
                                                             
                                                             1)), k - 1)
            
            else addprior <- rep(prior.param[i], k - 1)
            
          }
          
          new.prior <- c(new.prior, addprior)
          
        }
        
        prior.param <- new.prior
        
        x <- leaps.x <- mmm
        
      }
      
    }
    
    
    
    
    
    xx <- data.frame()
    
    xx <- dropcols(leaps.x, surv.t, cens, maxCol)
    
    var.names <- xx$var.names
    
    remaining <- xx$remaining.vars
    
    leaps.x <- xx$mm
    
    reduced <- xx$any.dropped
    
    dropped <- NULL
    
    if (reduced) 
      
      dropped <- xx$dropped
    
    nvar <- length(x[1, ])
    
    x <- x[, remaining, drop = FALSE]
    
    x <- data.frame(x)
    
    fac.levels <- fac.levels[remaining]
    
    output.names <- list()
    
    for (i in 1:length(var.names)) {
      
      if (is.factor(x[, i])) 
        
        output.names[[i]] <- levels(x[, i])
      
      else output.names[[i]] <- NA
      
    }
    
    xnames <- names(x)
    
    names(leaps.x) <- var.names
    
    x.coxph <- data.frame(surv.t = surv.t, cens = cens, leaps.x)
    
    cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                     
                     method = "breslow")
    
    x.coxph.fac <- data.frame(surv.t = surv.t, cens = cens, x)
    
    cox.assign <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph.fac, 
                        
                        method = "breslow")$assign
    
    for (tempi in length(cox.assign):1) {
      
      cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
      
      names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
      
    }
    
    cox.assign[[1]] <- numeric(0)
    
    names(cox.assign)[1] <- "(Intercept)"
    
    n <- sum(cens)
    
    if (ncol(leaps.x) >= 3) {
      
      coef <- cox.out$coef
      
      info <- solve(cox.out$var)
      
      a <- leaps.bs(info, coef, names.arg = xnames, nbest = nbest)
      
      a$r2 <- c(0, a$r2)
      
      a$size <- c(0, a$size)
      
      a$label <- c("NULL", a$label)
      
      a$which <- rbind(rep(FALSE, ncol(x)), a$which)
      
      nmod <- length(a$size)
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(a$which * prior.mat + (!a$which) * (1 - 
                                                           
                                                           prior.mat), 1, prod)
      
      bIb <- as.numeric(coef %*% info %*% coef)
      
      lrt <- bIb - (a$r2 * bIb)
      
      aic <- lrt + (a$size) * log(n) - 2 * log(prior)
      
      occam <- aic - min(aic) < 2 * OR.fix * log(OR)
      
      size <- a$size[occam]
      
      label <- a$label[occam]
      
      which <- a$which[occam, , drop = FALSE]
      
      aic <- aic[occam]
      
      prior <- prior[occam]
      
    }
    
    else {
      
      nmod <- switch(ncol(x), 2, 4)
      
      aic <- label <- rep(0, nmod)
      
      which <- matrix(c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE), nmod, nmod/2)
      
      size <- c(0, 1, 1, 2)[1:nmod]
      
      sep <- ","
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(which * prior.mat + (!which) * (1 - prior.mat), 
                     
                     1, prod)
      
      for (k in 1:nmod) {
        
        if (k == 1) 
          
          label[k] <- "NULL"
        
        else label[k] <- paste(names(leaps.x)[which[k, ]], 
                               
                               collapse = sep)
        
      }
      
    }
    
    model.fits <- as.list(rep(0, length(label)))
    
    loglik <- rep(0, length(label))
    
    size <- rep(0, length(label))
    
    loglik.null <- cox.out$loglik[1]
    
    for (k in (1:length(label))) {
      
      if (sum(which[k, ]) != 0) {
        
        x.coxph <- data.frame(x[, which[k, ], drop = FALSE], 
                              
                              surv.t = surv.t, cens = cens)
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         iter.max = 30, method = "breslow")
        
        loglik[k] <- cox.out$loglik[2]
        
        size[k] <- length(cox.out$coef)
        
        model.fits[[k]] <- matrix(rep(0, 2 * length(cox.out$coef)), 
                                  
                                  ncol = 2)
        
        model.fits[[k]][, 1] <- cox.out$coef
        
        model.fits[[k]][, 2] <- sqrt(diag(cox.out$var))
        
      }
      
      else {
        
        loglik[k] <- loglik.null
        
      }
      
    }
    
    aic <- 2 * size - 2 * (loglik - loglik.null) - 2 * log(prior)
    
    occam <- aic - min(aic) < 2 * log(OR)
    
    size <- size[occam]
    
    label <- label[occam]
    
    which <- which[occam, , drop = FALSE]
    
    aic <- aic[occam]
    
    prior <- prior[occam]
    
    model.fits <- model.fits[occam]
    
    postprob <- (exp(-0.5 * (aic - min(aic))))/sum(exp(-0.5 * 
                                                         
                                                         (aic - min(aic))))
    
    order.aic <- order(aic, size, label)
    
    size <- size[order.aic]
    
    label <- label[order.aic]
    
    which <- which[order.aic, , drop = FALSE]
    
    aic <- aic[order.aic]
    
    prior <- prior[order.aic]
    
    postprob <- postprob[order.aic]
    
    model.fits <- model.fits[order.aic]
    
    nmod <- length(size)
    
    if (strict & (nmod != 1)) {
      
      occam <- rep(TRUE, nmod)
      
      for (k in (2:nmod)) {
        
        for (j in (1:(k - 1))) {
          
          which.diff <- which[k, ] - which[j, ]
          
          if (all(which.diff >= 0)) 
            
            occam[k] <- FALSE
          
        }
        
      }
      
      size <- size[occam]
      
      label <- label[occam]
      
      which <- which[occam, , drop = FALSE]
      
      aic <- aic[occam]
      
      prior <- prior[occam]
      
      postprob <- postprob[occam]/(sum(postprob[occam]))
      
      model.fits <- model.fits[occam]
      
    }
    
    probne0 <- round(100 * t(which) %*% as.matrix(postprob), 
                     
                     1)
    
    nmod <- length(aic)
    
    nvar <- max(unlist(cox.assign))
    
    Ebi <- rep(0, nvar)
    
    SDbi <- rep(0, nvar)
    
    EbiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    sebiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    for (i in (1:ncol(x))) {
      
      whereisit <- cox.assign[[i + 1]]
      
      if (any(which[, i])) 
        
        for (k in (1:nmod)) if (which[k, i] == TRUE) {
          
          spot <- sum(which[k, (1:i)])
          
          posMk <- (c(0, cumsum(fac.levels[which[k, ]])) + 
                      
                      1)[spot]
          
          posMk <- posMk:(posMk + fac.levels[i] - 1)
          
          EbiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                 
                                                 1]
          
          sebiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                  
                                                  2]
          
        }
      
    }
    
    Ebi <- postprob %*% EbiMk
    
    Ebimat <- matrix(rep(Ebi, nmod), nrow = nmod, byrow = TRUE)
    
    SDbi <- sqrt(postprob %*% (sebiMk^2) + postprob %*% ((EbiMk - 
                                                            
                                                            Ebimat)^2))
    
    CSDbi <- rep(0, nvar)
    
    CEbi <- CSDbi
    
    for (i in (1:ncol(x))) {
      
      sel <- which[, i]
      
      if (sum(sel) > 0) {
        
        cpp <- rbind(postprob[sel]/sum(postprob[sel]))
        
        CEbi[cox.assign[[i + 1]]] <- as.numeric(cpp %*% EbiMk[sel, 
                                                              
                                                              cox.assign[[i + 1]]])
        
        CSDbi[cox.assign[[i + 1]]] <- sqrt(cpp %*% (sebiMk[sel, 
                                                           
                                                           cox.assign[[i + 1]]]^2) + cpp %*% ((EbiMk[sel, 
                                                                                                     
                                                                                                     cox.assign[[i + 1]]] - CEbi[cox.assign[[i + 1]]])^2))
        
      }
      
    }
    
    names(output.names) <- var.names
    
    result <- list(postprob = postprob, label = label, size = size, 
                   
                   aic = aic, prior.param = prior.param, prior.model.weights = prior/prior.weight.denom, 
                   
                   which = which, probne0 = c(probne0), postmean = as.vector(Ebi), 
                   
                   postsd = as.vector(SDbi), condpostmean = CEbi, condpostsd = CSDbi, 
                   
                   mle = EbiMk, se = sebiMk, namesx = var.names, reduced = reduced, 
                   
                   dropped = dropped, call = cl, n.models = length(postprob), 
                   
                   n.vars = length(probne0), nests = length(Ebi), output.names = output.names, 
                   
                   assign = cox.assign, factor.type = factor.type)
    
    class(result) <- "aic.surv"
    
    result
    
  }

"aic.surv.formula" <-
  
  function (f, data, strict = FALSE, OR = 20, maxCol = 30, prior.param = c(rep(0.5, 
                                                                               
                                                                               ncol(x))), OR.fix = 2, nbest = 150, factor.type = TRUE, factor.prior.adjust = FALSE, 
            
            call = NULL, ...) 
    
  {
    
    cl <- match.call()
    
    tms <- terms(f, data = data)
    
    fmatrix <- attr(tms, "factors")
    
    tms.order <- attr(tms, "order")
    
    tms.labels <- attr(tms, "term.labels")
    
    mm <- model.matrix(tms, data = data)
    
    assn <- attr(mm, "assign")
    
    nterms <- max(assn)
    
    datalist <- eval(attr(tms, "variables"), envir = data)
    
    nvar <- nrow(fmatrix) - 1
    
    isvarfac <- rep(NA, times = nvar)
    
    for (i in 1:nvar) isvarfac[i] <- is.factor(datalist[[i + 
                                                           
                                                           1]])
    
    istermfac <- rep(NA, times = nterms)
    
    for (i in 1:nterms) {
      
      cterms <- fmatrix[-1, i] == 1
      
      istermfac[i] <- sum(isvarfac[cterms] == FALSE) == 0
      
    }
    
    surv.t.name <- all.vars(f)[1]
    
    cens.name <- all.vars(f)[2]
    
    moddata <- data.frame(rep(NA, times = dim(mm)[1]))
    
    cnames <- NULL
    
    for (i in 1:nterms) {
      
      if (istermfac[i]) {
        
        if (tms.order[i] == 1) {
          
          moddata <- cbind(moddata, datalist[[i + 1]])
          
          cnames <- c(cnames, tms.labels[i])
          
        }
        
        else {
          
          sel <- assn == i
          
          nlev <- sum(sel)
          
          newfac.index <- (mm[, sel] %*% cbind(1:nlev)) + 
            
            1
          
          facnames <- c("ref", colnames(mm)[sel])
          
          newfac <- facnames[newfac.index]
          
          newfac <- factor(newfac)
          
          moddata <- cbind(moddata, newfac)
          
          cnames <- c(cnames, paste(tms.labels[i], "..", 
                                    
                                    sep = ""))
          
        }
        
      }
      
      else {
        
        sel <- assn == i
        
        moddata <- cbind(moddata, mm[, sel])
        
        cnames <- c(cnames, colnames(mm)[sel])
        
      }
      
    }
    
    moddata <- moddata[, -1]
    
    cnames <- gsub(":", ".", cnames)
    
    moddata <- cbind(moddata, datalist[[1]][, 1], datalist[[1]][, 
                                                                
                                                                2])
    
    colnames(moddata) <- c(cnames, surv.t.name, cens.name)
    
    nv <- ncol(moddata) - 2
    
    surv.t <- moddata[, nv + 1]
    
    cens <- moddata[, nv + 2]
    
    x <- moddata[, 1:nv]
    
    aic.surv(x, surv.t, cens, strict = FALSE, OR = OR, maxCol = maxCol, 
             
             prior.param = prior.param, OR.fix = OR.fix, nbest = nbest, 
             
             factor.type = factor.type, factor.prior.adjust = factor.prior.adjust, 
             
             call = cl)
    
  }

"aic.surv.matrix" <-
  
  function (x, surv.t, cens, strict = FALSE, OR = 20, maxCol = 30, 
            
            prior.param = c(rep(0.5, ncol(x))), OR.fix = 2, nbest = 150, 
            
            factor.type = TRUE, factor.prior.adjust = FALSE, call = NULL, ...) 
    
  {
    
    leaps.bs <- function(info, coef, names.arg, nbest = nbest) {
      
      names.arg <- names.arg
      
      if (is.null(names.arg)) 
        
        names.arg <- c(as.character(1:9), LETTERS, letters)[1:ncol(info)]
      
      if (length(names.arg) < ncol(info)) 
        
        stop("Too few names")
      
      bIb <- coef %*% info %*% coef
      
      kx <- ncol(info)
      
      maxreg <- nbest * kx
      
      if (kx < 3) 
        
        stop("Too few independent variables")
      
      imeth <- 1
      
      df <- kx + 1
      
      Ib <- info %*% coef
      
      rr <- cbind(info, Ib)
      
      rr <- rbind(rr, c(Ib, bIb))
      
      it <- 0
      
      n.cols <- kx + 1
      
      nv <- kx + 1
      
      nf <- 0
      
      no <- 1e+05
      
      ib <- 1
      
      mb <- nbest
      
      nd <- n.cols
      
      nc <- 4 * n.cols
      
      rt <- matrix(rep(0, times = nd * nc), ncol = nc)
      
      rt[, 1:n.cols] <- rr
      
      iw <- c(1:(kx + 1), rep(0, times = 4 * nd))
      
      nw <- length(iw)
      
      rw <- rep(0, times = 2 * mb * kx + 7 * nd)
      
      nr <- length(rw)
      
      t1 <- 2
      
      s2 <- -1
      
      ne <- 0
      
      iv <- 0
      
      nret <- mb * kx
      
      Subss <- rep(0, times = nret)
      
      RSS <- Subss
      
      ans <- .Fortran("fwleaps", as.integer(nv), as.integer(it), 
                      
                      as.integer(kx), as.integer(nf), as.integer(no), as.integer(1), 
                      
                      as.double(2), as.integer(mb), as.double(rt), as.integer(nd), 
                      
                      as.integer(nc), as.integer(iw), as.integer(nw), as.double(rw), 
                      
                      as.integer(nr), as.double(t1), as.double(s2), as.integer(ne), 
                      
                      as.integer(iv), as.double(Subss), as.double(RSS), 
                      
                      as.integer(nret), PACKAGE = "BMA")
      
      regid <- ans[[21]]/2
      
      r2 <- ans[[20]]
      
      nreg <- sum(regid > 0)
      
      regid <- regid[1:nreg]
      
      r2 <- r2[1:nreg]
      
      which <- matrix(TRUE, nreg, kx)
      
      z <- regid
      
      which <- matrix(as.logical((rep.int(z, kx)%/%rep.int(2^((kx - 
                                                                 
                                                                 1):0), rep.int(length(z), kx)))%%2), byrow = FALSE, ncol = kx)
      
      size <- which %*% rep(1, kx)
      
      label <- character(nreg)
      
      sep <- if (all(nchar(names.arg) == 1)) 
        
        ""
      
      else ","
      
      for (i in 1:nreg) label[i] <- paste(names.arg[which[i, 
                                                          
                                                          ]], collapse = sep)
      
      ans <- list(r2 = r2, size = size, label = label, which = which)
      
      return(ans)
      
    }
    
    factor.names <- function(x) {
      
      out <- list()
      
      for (i in 1:ncol(x)) if (is.factor(x[, i])) 
        
        out[[i]] <- levels(x[, i])
      
      else out <- c(out, list(NULL))
      
      attributes(out)$names <- names(x)
      
      return(out)
      
    }
    
    dropcols <- function(x, surv.t, cens, maxCols = 30) {
      
      vnames <- attributes(x)$names
      
      nvar <- length(vnames)
      
      isfac <- rep(FALSE, times = nvar)
      
      for (i in 1:nvar) isfac[i] <- is.factor(x[, i])
      
      nlevels <- rep(NA, times = nvar)
      
      for (i in 1:nvar) if (isfac[i]) 
        
        nlevels[i] <- length(levels(x[, i]))
      
      any.dropped <- FALSE
      
      mm <- model.matrix(terms.formula(~., data = x), data = x)
      
      designx <- attributes(mm)$assign
      
      n.designx <- length(designx)
      
      designx.levels <- rep(1, times = n.designx)
      
      for (i in 2:n.designx) if (isfac[designx[i]]) 
        
        designx.levels[i] <- sum(designx[1:i] == designx[i]) + 
        
        1
      
      x.coxph <- data.frame(mm[, -1], surv.t = surv.t, cens = cens)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow", iter.max = 30)
      
      while (length(cox.out$coefficients) > maxCol) {
        
        any.dropped <- TRUE
        
        dropcox <- drop1(cox.out, test = "Chisq")
        
        dropped <- which.max(dropcox$"Pr(Chi)"[-1]) + 1
        
        x.coxph <- x.coxph[, -(dropped - 1)]
        
        designx.levels <- designx.levels[-dropped]
        
        designx <- designx[-dropped]
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         method = "breslow", iter.max = 30)
        
      }
      
      remaining.vars <- unique(designx[-1])
      
      new.nvar <- length(remaining.vars)
      
      dropped.vars <- vnames[-remaining.vars]
      
      dropped.levels <- NULL
      
      ncol.cox <- ncol(x.coxph) - 2
      
      x.coxph <- x.coxph[-((ncol.cox + 1):(ncol.cox + 2))]
      
      xx <- data.frame(matrix(rep(NA, times = new.nvar * nrow(x.coxph)), 
                              
                              ncol = new.nvar))
      
      new.names = rep(NA, times = new.nvar)
      
      for (i in 1:new.nvar) {
        
        cvar <- remaining.vars[i]
        
        lvls <- designx.levels[cvar == designx]
        
        if (isfac[cvar]) {
          
          if (length(lvls) != length(levels(x[, cvar]))) {
            
            newvar <- (as.matrix(x.coxph[, cvar == designx[-1]]) %*% 
                         
                         cbind(lvls - 1)) + 1
            
            xx[, i] <- factor(levels(x[, cvar])[newvar])
            
            new.names[i] <- vnames[cvar]
            
            removed.levels <- levels(x[, cvar])[-c(1, lvls)]
            
            dropped.levels <- c(dropped.levels, paste(vnames[cvar], 
                                                      
                                                      "_", removed.levels, sep = ""))
            
          }
          
          else {
            
            xx[, i] <- factor(x[, cvar])
            
            new.names[i] <- vnames[cvar]
            
          }
          
        }
        
        else {
          
          xx[, i] <- x[, cvar]
          
          new.names[i] <- vnames[cvar]
          
        }
        
      }
      
      dropped <- c(dropped.vars, dropped.levels)
      
      return(list(mm = xx, any.dropped = any.dropped, dropped = dropped, 
                  
                  var.names = new.names, remaining.vars = remaining.vars))
      
    }
    
    if (is.null(call)) 
      
      cl <- match.call()
    
    else cl <- call
    
    options(contrasts = c("contr.treatment", "contr.treatment"))
    
    x <- data.frame(x)
    
    prior.weight.denom <- 0.5^ncol(x)
    
    x.omit <- na.omit(x)
    
    used <- match(row.names(x), row.names(x.omit))
    
    omitted <- seq(nrow(x))[is.na(used)]
    
    if (length(omitted) > 0) {
      
      x <- x.omit
      
      surv.t <- surv.t[-omitted]
      
      cens <- cens[-omitted]
      
      warning(paste("There were ", length(omitted), "records deleted due to NA's"))
      
    }
    
    leaps.x <- x
    
    fac.levels <- rep(1, times = ncol(x))
    
    fn <- factor.names(x)
    
    factors <- !all(unlist(lapply(fn, is.null)))
    
    if (factors) {
      
      cdf <- cbind.data.frame(y = surv.t, x)
      
      mm <- model.matrix(formula(cdf), data = cdf)[, -1, drop = FALSE]
      
      mmm <- data.frame(matrix(mm, nrow = nrow(mm), byrow = FALSE))
      
      names(mmm) <- dimnames(mm)[[2]]
      
      output.names <- names(mmm)
      
      x.coxph <- data.frame(surv.t = surv.t, cens = cens, x)
      
      cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                       
                       method = "breslow")
      
      cox.assign <- cox.out$assign
      
      for (tempi in length(cox.assign):1) {
        
        cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
        
        names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
        
      }
      
      cox.assign[[1]] <- numeric(0)
      
      names(cox.assign)[1] <- "(Intercept)"
      
      fac.levels <- unlist(lapply(cox.assign, length)[-1])
      
      if (factor.type) {
        
        for (i in 1:length(names(x))) {
          
          if (!is.null(fn[[i]])) {
            
            nx <- names(x)[i]
            
            coefs <- cox.out$coef[cox.out$assign[[i]]]
            
            old.vals <- x[, i]
            
            new.vals <- c(0, coefs)
            
            new.vec <- as.vector(new.vals[match(old.vals, 
                                                
                                                fn[[i]])])
            
            leaps.x[, nx] <- new.vec
            
          }
          
        }
        
      }
      
      else {
        
        fac.levels <- rep(1, times = ncol(mmm))
        
        new.prior <- NULL
        
        for (i in 1:length(names(x))) {
          
          addprior <- prior.param[i]
          
          if (!is.null(fn[[i]])) {
            
            k <- length(fn[[i]])
            
            if (factor.prior.adjust) 
              
              addprior <- rep(1 - (1 - prior.param[i])^(1/(k - 
                                                             
                                                             1)), k - 1)
            
            else addprior <- rep(prior.param[i], k - 1)
            
          }
          
          new.prior <- c(new.prior, addprior)
          
        }
        
        prior.param <- new.prior
        
        x <- leaps.x <- mmm
        
      }
      
    }
    
    xx <- data.frame()
    
    xx <- dropcols(leaps.x, surv.t, cens, maxCol)
    
    var.names <- xx$var.names
    
    remaining <- xx$remaining.vars
    
    leaps.x <- xx$mm
    
    reduced <- xx$any.dropped
    
    dropped <- NULL
    
    if (reduced) 
      
      dropped <- xx$dropped
    
    nvar <- length(x[1, ])
    
    x <- x[, remaining, drop = FALSE]
    
    x <- data.frame(x)
    
    fac.levels <- fac.levels[remaining]
    
    output.names <- list()
    
    for (i in 1:length(var.names)) {
      
      if (is.factor(x[, i])) 
        
        output.names[[i]] <- levels(x[, i])
      
      else output.names[[i]] <- NA
      
    }
    
    xnames <- names(x)
    
    names(leaps.x) <- var.names
    
    x.coxph <- data.frame(surv.t = surv.t, cens = cens, leaps.x)
    
    cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                     
                     method = "breslow")
    
    x.coxph.fac <- data.frame(surv.t = surv.t, cens = cens, x)
    
    cox.assign <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph.fac, 
                        
                        method = "breslow")$assign
    
    for (tempi in length(cox.assign):1) {
      
      cox.assign[[tempi + 1]] <- cox.assign[[tempi]]
      
      names(cox.assign)[tempi + 1] <- names(cox.assign)[tempi]
      
    }
    
    cox.assign[[1]] <- numeric(0)
    
    names(cox.assign)[1] <- "(Intercept)"
    
    n <- sum(cens)
    
    if (ncol(leaps.x) >= 3) {
      
      coef <- cox.out$coef
      
      info <- solve(cox.out$var)
      
      a <- leaps.bs(info, coef, names.arg = xnames, nbest = nbest)
      
      a$r2 <- c(0, a$r2)
      
      a$size <- c(0, a$size)
      
      a$label <- c("NULL", a$label)
      
      a$which <- rbind(rep(FALSE, ncol(x)), a$which)
      
      nmod <- length(a$size)
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(a$which * prior.mat + (!a$which) * (1 - 
                                                           
                                                           prior.mat), 1, prod)
      
      bIb <- as.numeric(coef %*% info %*% coef)
      
      lrt <- bIb - (a$r2 * bIb)
      
      aic <- lrt + (a$size) * log(n) - 2 * log(prior)
      
      occam <- aic - min(aic) < 2 * OR.fix * log(OR)
      
      size <- a$size[occam]
      
      label <- a$label[occam]
      
      which <- a$which[occam, , drop = FALSE]
      
      aic <- aic[occam]
      
      prior <- prior[occam]
      
    }
    
    else {
      
      nmod <- switch(ncol(x), 2, 4)
      
      aic <- label <- rep(0, nmod)
      
      which <- matrix(c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE), nmod, nmod/2)
      
      size <- c(0, 1, 1, 2)[1:nmod]
      
      sep <- ","
      
      prior.mat <- matrix(rep(prior.param, nmod), nmod, ncol(x), 
                          
                          byrow = TRUE)
      
      prior <- apply(which * prior.mat + (!which) * (1 - prior.mat), 
                     
                     1, prod)
      
      for (k in 1:nmod) {
        
        if (k == 1) 
          
          label[k] <- "NULL"
        
        else label[k] <- paste(names(leaps.x)[which[k, ]], 
                               
                               collapse = sep)
        
      }
      
    }
    
    model.fits <- as.list(rep(0, length(label)))
    
    loglik <- rep(0, length(label))
    
    size <- rep(0, length(label))
    
    loglik.null <- cox.out$loglik[1]
    
    for (k in (1:length(label))) {
      
      if (sum(which[k, ]) != 0) {
        
        x.coxph <- data.frame(x[, which[k, ], drop = FALSE], 
                              
                              surv.t = surv.t, cens = cens)
        
        cox.out <- coxph(Surv(surv.t-1, surv.t, cens) ~ ., data = x.coxph, 
                         
                         iter.max = 30, method = "breslow")
        
        loglik[k] <- cox.out$loglik[2]
        
        size[k] <- length(cox.out$coef)
        
        model.fits[[k]] <- matrix(rep(0, 2 * length(cox.out$coef)), 
                                  
                                  ncol = 2)
        
        model.fits[[k]][, 1] <- cox.out$coef
        
        model.fits[[k]][, 2] <- sqrt(diag(cox.out$var))
        
        aic[k] <- AIC(cox.out)
        
      }
      
      else {
        
        loglik[k] <- loglik.null
        
      }
      
    }
    
    occam <- aic - min(aic) < 2 * log(OR)
    
    size <- size[occam]
    
    label <- label[occam]
    
    which <- which[occam, , drop = FALSE]
    
    aic <- aic[occam]
    
    prior <- prior[occam]
    
    model.fits <- model.fits[occam]
    
    postprob <- (exp(-0.5 * (aic - min(aic))))/sum(exp(-0.5 * 
                                                         
                                                         (aic - min(aic))))
    
    order.aic <- order(aic, size, label)
    
    size <- size[order.aic]
    
    label <- label[order.aic]
    
    which <- which[order.aic, , drop = FALSE]
    
    aic <- aic[order.aic]
    
    prior <- prior[order.aic]
    
    postprob <- postprob[order.aic]
    
    model.fits <- model.fits[order.aic]
    
    nmod <- length(size)
    
    if (strict & (nmod != 1)) {
      
      occam <- rep(TRUE, nmod)
      
      for (k in (2:nmod)) {
        
        for (j in (1:(k - 1))) {
          
          which.diff <- which[k, ] - which[j, ]
          
          if (all(which.diff >= 0)) 
            
            occam[k] <- FALSE
          
        }
        
      }
      
      size <- size[occam]
      
      label <- label[occam]
      
      which <- which[occam, , drop = FALSE]
      
      aic <- aic[occam]
      
      prior <- prior[occam]
      
      postprob <- postprob[occam]/(sum(postprob[occam]))
      
      model.fits <- model.fits[occam]
      
    }
    
    probne0 <- round(100 * t(which) %*% as.matrix(postprob), 
                     
                     1)
    
    nmod <- length(aic)
    
    nvar <- max(unlist(cox.assign))
    
    Ebi <- rep(0, nvar)
    
    SDbi <- rep(0, nvar)
    
    EbiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    sebiMk <- matrix(rep(0, nmod * nvar), nrow = nmod)
    
    for (i in (1:ncol(x))) {
      
      whereisit <- cox.assign[[i + 1]]
      
      if (any(which[, i])) 
        
        for (k in (1:nmod)) if (which[k, i] == TRUE) {
          
          spot <- sum(which[k, (1:i)])
          
          posMk <- (c(0, cumsum(fac.levels[which[k, ]])) + 
                      
                      1)[spot]
          
          posMk <- posMk:(posMk + fac.levels[i] - 1)
          
          EbiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                 
                                                 1]
          
          sebiMk[k, whereisit] <- model.fits[[k]][posMk, 
                                                  
                                                  2]
          
        }
      
    }
    
    Ebi <- postprob %*% EbiMk
    
    Ebimat <- matrix(rep(Ebi, nmod), nrow = nmod, byrow = TRUE)
    
    SDbi <- sqrt(postprob %*% (sebiMk^2) + postprob %*% ((EbiMk - 
                                                            
                                                            Ebimat)^2))
    
    CSDbi <- rep(0, nvar)
    
    CEbi <- CSDbi
    
    for (i in (1:ncol(x))) {
      
      sel <- which[, i]
      
      if (sum(sel) > 0) {
        
        cpp <- rbind(postprob[sel]/sum(postprob[sel]))
        
        CEbi[cox.assign[[i + 1]]] <- as.numeric(cpp %*% EbiMk[sel, 
                                                              
                                                              cox.assign[[i + 1]]])
        
        CSDbi[cox.assign[[i + 1]]] <- sqrt(cpp %*% (sebiMk[sel, 
                                                           
                                                           cox.assign[[i + 1]]]^2) + cpp %*% ((EbiMk[sel, 
                                                                                                     
                                                                                                     cox.assign[[i + 1]]] - CEbi[cox.assign[[i + 1]]])^2))
        
      }
      
    }
    
    names(output.names) <- var.names
    
    result <- list(postprob = postprob, label = label, size = size, 
                   
                   aic = aic, prior.param = prior.param, prior.model.weights = prior/prior.weight.denom, 
                   
                   which = which, probne0 = c(probne0), postmean = as.vector(Ebi), 
                   
                   postsd = as.vector(SDbi), condpostmean = CEbi, condpostsd = CSDbi, 
                   
                   mle = EbiMk, se = sebiMk, namesx = var.names, reduced = reduced, 
                   
                   dropped = dropped, call = cl, n.models = length(postprob), 
                   
                   n.vars = length(probne0), nests = length(Ebi), output.names = output.names, 
                   
                   assign = cox.assign, factor.type = factor.type)
    
    class(result) <- "aic.surv"
    
    result
    
  }

"summary.aic.surv" <-
  
  function (object, n.models = 5, digits = max(3, getOption("digits") - 
                                                 
                                                 3), conditional = FALSE, display.dropped = FALSE, ...) 
    
  {
    
    x<- object
    
    cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
    
    if (display.dropped & x$reduced) {
      
      cat("\nThe following variables were dropped prior to averaging:\n")
      
      cat(x$dropped)
      
      cat("\n")
      
    }
    
    n.models <- min(n.models, x$n.models)
    
    sel <- 1:n.models
    
    cat("\n ", length(x$postprob), " models were selected")
    
    cat("\n Best ", n.models, " models (cumulative posterior probability = ", 
        
        round(sum(x$postprob[sel]), digits), "): \n\n")
    
    nms <- length(x$namesx)
    
    ncx <- length(unlist(x$assign)[-1])
    
    nvar <- rep(0, times = n.models)
    
    for (i in 1:nms) nvar <- nvar + as.numeric(as.vector(rbind(rep(1, 
                                                                   
                                                                   length(x$assign[[i + 1]]))) %*% (t(x$mle[sel, x$assign[[i + 
                                                                                                                             
                                                                                                                             1]], drop = FALSE] != 0)) > 0))
    
    modelposts <- format(round(x$postprob[sel], 3), digits = 3)
    
    coeffs <- t(x$mle[sel, , drop = FALSE])
    
    cfaic <- rbind(x$aic[sel], coeffs)
    
    cfaicf <- format(cfaic, digits = digits)
    
    coeffsf <- cfaicf[-1, , drop = FALSE]
    
    aic <- cfaicf[1, , drop = FALSE]
    
    postmeans <- format(x$postmean, digits = digits)
    
    postsds <- format(x$postsd, digits = digits)
    
    if (conditional) {
      
      cpostmeans <- format(x$condpostmean, digits = digits)
      
      cpostsds <- format(x$condpostsd, digits = digits)
      
    }
    
    varposts <- format(round(x$probne0, 1), digits = 3)
    
    strlength <- nchar(coeffsf[1, 1])
    
    decpos <- nchar(unlist(strsplit(coeffsf[1, 1], "\\."))[1])
    
    offset <- paste(rep(" ", times = decpos - 1), sep = "", collapse = "")
    
    offset2 <- paste(rep(" ", times = decpos + 1), sep = "", 
                     
                     collapse = "")
    
    modelposts <- paste(offset, modelposts, sep = "")
    
    nvar <- paste(offset2, nvar, sep = "")
    
    dotoffset <- round(max(nchar(coeffsf))/2)
    
    zerocoefstring <- paste(paste(rep(" ", times = dotoffset), 
                                  
                                  collapse = "", sep = ""), ".", sep = "")
    
    coeffsf[coeffs == 0] <- zerocoefstring
    
    avp <- NULL
    
    for (i in 1:nms) {
      
      avp <- rbind(avp, varposts[i])
      
      if (!is.na(x$output.names[[i]][1])) 
        
        avp <- rbind(avp, cbind(rep("", times = length(x$assign[[i + 
                                                                   
                                                                   1]]))))
      
    }
    
    top <- cbind(postmeans, postsds)
    
    if (conditional) 
      
      top <- cbind(top, cpostmeans, cpostsds)
    
    top <- cbind(top, coeffsf)
    
    atop <- NULL
    
    for (i in 1:nms) {
      
      if (!is.na(x$output.names[[i]][1])) 
        
        atop <- rbind(atop, rbind(rep("", times = ncol(top))))
      
      atop <- rbind(atop, top[x$assign[[i + 1]], ])
      
    }
    
    top <- cbind(avp, atop)
    
    linesep <- rep("", times = ncol(top))
    
    offset <- c("", "", "")
    
    if (conditional) 
      
      offset <- c(offset, c("", ""))
    
    bottom <- rbind(c(offset, nvar), c(offset, aic), c(offset, 
                                                       
                                                       modelposts))
    
    out <- rbind(top, linesep, bottom)
    
    vnames <- NULL
    
    for (i in 1:nms) {
      
      vnames <- c(vnames, names(x$output.names[i]))
      
      blnk <- paste(rep(" ", times = nchar(names(x$output.names[i]))), 
                    
                    collapse = "")
      
      if (!is.na(x$output.names[i][1])) 
        
        vnames <- c(vnames, paste(blnk, unlist(x$output.names[i])[-1], 
                                  
                                  sep = "."))
      
    }
    
    row.names(out) <- c(vnames, "", "nVar", "AIC", "post prob")
    
    colnms <- c("p!=0", " EV", "SD")
    
    if (conditional) 
      
      colnms <- c(colnms, "cond EV", "cond SD")
    
    colnms <- c(colnms, paste("model ", 1:n.models, sep = ""))
    
    dimnames(out)[[2]] <- colnms
    
    print.default(out, print.gap = 2, quote = FALSE, ...)
    
  }

"plot.bic.surv" <-
  
  function (x, e = 1e-04, mfrow = NULL, include = 1:length(x$namesx), ...) 
    
  {
    
    plotvar <- function(prob0, mixprobs, means, sds, Emean, Esd, 
                        
                        name, e = 1e-04, nsteps = 500, ...) {
      
      if (prob0 == 1) {
        
        xlower <- -0
        
        xupper <- 0
        
        xmax <- 1
        
      }
      
      else {
        
        qmin <- qnorm(e/2, Emean, Esd)
        
        qmax = qnorm(1 - e/2, Emean, Esd)
        
        xlower <- min(qmin, 0)
        
        xupper <- max(0, qmax)
        
      }
      
      xx <- seq(xlower, xupper, length.out = nsteps)
      
      yy <- rep(0, times = length(xx))
      
      maxyy <- 1
      
      if (prob0 < 1) {
        
        for (j in 1:length(means)) yy <- yy + mixprobs[j] * 
            
            dnorm(xx, means[j], sds[j])
        
        maxyy <- max(yy)
        
      }
      
      ymax <- max(prob0, 1 - prob0)
      
      plot(c(xlower, xupper), c(0, ymax), type = "n", xlab = "", 
           
           ylab = "", main = name, ...)
      
      lines(c(0, 0), c(0, prob0), lty = 1, lwd = 3, ...)
      
      lines(xx, (1 - prob0) * yy/maxyy, lty = 1, lwd = 1, ...)
      
    }
    
    vars <- unlist(x$assign[include + 1])
    
    nvar <- length(vars)
    
    probs <- NULL
    
    for (i in include) probs <- c(probs, rep(x$probne0[i], times = length(x$assign[[i + 
                                                                                      
                                                                                      1]])))
    
    nms <- NULL
    
    for (i in include) {
      
      if (is.na(x$output.names[i][1])) 
        
        nms <- c(nms, names(x$output.names[i]))
      
      else nms <- c(nms, paste(names(x$output.names[i]), unlist(x$output.names[i])[-1], 
                               
                               sep = "."))
      
    }
    
    wwhich <- NULL
    
    for (i in include) {
      
      wwhich <- cbind(wwhich, matrix(rep(x$which[, i], times = length(x$assign[[i + 
                                                                                  
                                                                                  1]])), ncol = length(x$assign[[i + 1]])))
      
    }
    
    if (!is.null(mfrow)) {
      
      lo <- mfrow
      
      losize <- lo[1] * lo[2]
      
    }
    
    else {
      
      layoutsizes <- c(1, 4, 9)
      
      layouts <- rbind(c(1, 1), c(2, 2), c(3, 3))
      
      layout <- max((1:length(layoutsizes))[layoutsizes <= 
                                              
                                              nvar])
      
      losize <- layoutsizes[layout]
      
      lo <- layouts[layout, ]
      
    }
    
    
    
    keep.mfrow = par()$mfrow
    
    par(mfrow = c(1, 1))
    
    par(ask = FALSE)
    
    par(mfrow = lo)
    
    ngroups <- ceiling(nvar/losize)
    
    for (k in 1:ngroups) {
      
      for (ii in ((k - 1) * losize + 1):min(k * losize, nvar)) {
        
        i <- vars[ii]
        
        prob0 <- 1 - probs[i]/100
        
        sel <- wwhich[, i]
        
        mixprobs <- x$postprob[sel]/(1 - prob0)
        
        means <- x$mle[sel, i]
        
        sds <- x$se[sel, i]
        
        Emean <- x$condpostmean[i]
        
        Esd <- x$condpostsd[i]
        
        name <- nms[i]
        
        plotvar(prob0, mixprobs, means, sds, Emean, Esd, 
                
                name, e = e, ...)
        
      }
      
      par(ask = TRUE)
      
    }
    
    par(mfrow=keep.mfrow)
    
    par(ask = FALSE)
    
  }
